Quantum tunneling dynamics using hydrodynamic trajectories. 



Eric R. Bittner 

Department of Chemistry, University of Houston, Houston, TX 7720^ 
(Submitted to J. Chem. Phys.) 

In this paper we compute quantum trajectories arising from Bohm's causal description of quan- 
tum mechanics. Our computational methodology is based upon a finite-element moving least-squares 
method (MWLS) presented recently by Wyatt and co-workers (Lopreore and Wyatt, Phys. Rev. 
Lett. 82 , 5190 (1999)). This method treats the "particles" in the quantum Hamilton-Jacobi equa- 
tion as Lagrangian fluid elements which carry the phase, S, and density, p, required to reconstruct 
the quantum wavefunctions. Here, we compare results obtained via the MWLS procedure to ex- 
act results obtained either analytically or by numerical solution of the time dependent Schrodinger 
equation. Two systems are considered: firstly, dynamics in a harmonic well and secondly tunneling 
dynamics in a double well potential. In the case of tunneling in the double well potential, the 
quantum potential acts to lower the barrier separating the right and left hand sides of the well per- 
mitting trajectories to pass from one side to another. However, as probability density passes from 
one side to the other, the effective barrier begins to rise and eventually will segregate trajectories 
in one side from the other. We note that the MWLS trajectories exhibited long time stability in 
the purely harmonic cases. However, this stability was not evident in the barrier crossing dynamics. 
Comparisons to exact trajectories obtained via wave packet calculations indicate that the MWLS 
trajectories tend to underestimate the effects of constructive and destructive interference effects. 



I. INTRODUCTION 



Trajectory based constructions of quantum behavior are ubiquitous throughout quantum physics. In the semi- 
classical limit, quantum dynamics is approximated by classical equations of motion whereby the transition ampUtudes 
and wavefunctions are computed using the classical ftcAion connecting the initial and fWi^l states. Hydrodynamic 
constructions date back to early work by de BroglieErl3and Madelung' later by BohmJj'Q The so called quantum 
trajectories arising from this formalism have been the subject of a large number of ontological and philosophical 
papers seeking a causal interpretation of quantum mechanics. While a comprehensive overview of these works is far 
beyond the scope of this paper, Holland's book provides perhaps the most comprehensive technical overjaew of the 
approach and contains many examples of how the formalism can be applied in a wide variety of cases .□ In short, 
the quantum trajectories themselves are relatively easy to compute once one has obtained a wavefunction solution of 
Schrodinger equation, ij:. The velocity of a trajectory at a given point in space-time is computed as 

v[x,t)^——- (1) 
p{x,t} 

where j{x,t) is the quantum current, and p = In this "pilot wave" approach, ip acts as a guiding field for 

the trajectories. This approach is useful in that once the wave function has been obtained, is generally easy to 
compute the trajectories. While there are a number of papers which have computed quantum trajectories having 
obtained the wave function, relatively little work has been done in developiagi-Gpmputational methods based upon the 
description which does not rely upon first constructing the wave functionJHJ^Ej However, Wyatt and co-workers have 
recently described a mesh-less-SjUite-element method for integrating the de Broglie-Bohm equations using Lagrangian 
hydrodynamic trajectories.li3~E3 Similar methods are widely used in computational fluid dynamics (CFD) to, simulate 
fluid flow dynamics in porous media (such as an oil reservoir) and other systems with complex topologies.c3 Wyatt 's 
method represents the quantum density using a cloud of Lagrangian fluid points which themselves evolve according 
to the de Broglie-Bohm hydrodynamic equations of motion. The method features a moving weighted least-squares 
(MWLS) approach to compute the various derivatives and gradients required by the de Broglie-Bohm equations. 

In this work, we report on our implementation of the MWLS methodology and present an assessment of its difficulties 
and where improvements can be made. Previous applications of the approach have focused entirely upon reactive 
scattering type calculations. While this class of problems has its own set of associated difficulties, we have elected 
to focus upon the dynamics of systems trapped in various one dimensional potential wells. These problems, while 
perhaps too idealistic, allow one to compare trajectories and results obtained via a "hydrodynamic" calculation to 
those obtained via analytical or numerical solution of the time-dependent Schrodinger equation. Specifically, we first 
examine the dynamics of harmonic oscillators and then move onto a more challenging problem of tunneling in double 
well potential. Our results indicate that the methodology reproduces results obtained via more traditional grid based 
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approaches; however, we note that the long time stabihty of the methodology needs to be improved before it can be 
used as an alternative to wave-packet based calculation. 



II. QUANTUM EQUATIONS OF MOTION 

The de Broglie-Bohm equations are derived directly from the time dependent Schrodinger equation. The derivation 
is initialized by writing the wave function in polar form: 

V'(x,t) =i?(a;,i)e*^("'*)/'*. (2) 

When this substituted into the Schrodinger equation, the real and imaginary components can be collected into a 
continuity equation: 

dtp{x,t) + -V ■{pWS)^ (3) 
m 

and the quantum Hamilton-Jacobi equation: 

dtS+^-^^ + V + Q = (4) 

where Q is the so called quantum potential which is non-local and arises from the quantum kinetic energy operator 
in the Schrodinger equation. 

Q(.,t) = -|^^V^^/p (5) 

Taking the gradient of Eq.^, one obtains the equations of motion for the trajectories: 

Dtvix, t) - --V{y{x{t)) + Q{x{t))) (6) 
m 

where we define x{t) and v{x{t)) = V S/m as the trajectory and velocity of Lagrangian fluid elements. The notation 
Dtf denotes the "material" or "convective" derivative of the function / 

Dtf ^{dt + v-V)f (7) 

which gives the rate of change of f{{x{t)) as observed while moving with the particle along the trajectory, x{t) nil 
Hence in the Lagrangian picture of fluid mechanics, we "go with the flow" . The equation of continuity can also be 
written in terms of the convective derivative: 

Dtp+(y ■v)p^O (8) 

From this we can deduce a short time propagator for p: 

p{x,t + 6t)^er^-''^'p{x,t) (9) 

which is local in space and is dependent upon the divergence of the velocity field at the point x. To apply Eq. ^ 
and Eq. ^ in a computational scheme, we first discretize the system into an ensemble of Lagrangian fluid elements 
labeled by their position vectors, x. These elements carry information regarding the local density about that point 
and the phase is determined by integrating Eq, |^ along the path. From these two bits of information, the quantum 
wavefunction can be reconstructed by writingtS 

= \/p(a;i,0)exp (^^ L{x,,v^,T)dT + ^S{x,0)^ eiip (^-^ V • t;(r)dr j (10) 



where 



L{x,v,t) = ^mv^ - V -Q (11) 
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is the quantum Lagrangian and S{x, 0) represents any initial spatial dependent phase present in the initial wave 
function, such as S{x, 0) = Hkx. 

In order to compute expectation values we need to be able to integrate over the spatial domain spanned by either 
p or -tjj. If we require that each trajectory element carry a volume element, dx(t), then by normalization 
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1 = / dx{0)p{x, 0)^J2 P^iO)dx^{0) = P^it)dx^{t), (12) 

i=l 4=1 

where pi(t) — p(xi{t)) is the density carried by the ith trajectory. 

To compute dxi{t) we need the Jacobian which transforms the volume element dxi{0) to the volume element, dxi{t) 
some time later along trajectory Xi(t). 

Taking the material time derivative of J yields along each trajectory: 

AlogJ. ^V-w, (14) 
Thus, an amplitude element of the wavcfunction is propagated along a trajectory Xi(i)} as 

dxi{t)ipi{t) = dxi{0)y^pi{0)exp ^- J L{xi,Xi,T)dT + -S{xi,0)^ 



X exp y V • v,{T)dTj . (15) 

A typical propagation cycle consists of computing the quantum potential from the density at each fluid element, 
taking the gradient of Q and V to compute the acceleration of the fluid element and computing divergence of the 
velocity field to determine how the density in given fluid element changes over one discrete time interval. Thus, the 
trajectories, density, and velocity fields contain all of the essential information required to construct the full quantum 
mechanical wave function. 



III. MOVING LEAST SQUARES APPROXIMATION 



Since our goal is to be able to compute the quantum trajectories without computing the wavcfunction a priori, 
we need to be able to be able to discretize the system into a set of Lagrangian fluid elements, compute the various 
derivatives required to compute the equations of motion and flnally, advance the fluid elements to a new set of coordi- 
nates. This we accomplish through the use of an adaptive meshless^loud method described recently by WyattE3"E2. 
The method itself is described in detail by Liszka and co-workersE3 and is similar to methods used extensively in 
computational fluid dynamics. First, we review the finite element approximationEj (FEA) and its implementation in 
a movable weighted least-squares scheme (MWLS). 

In the moving least squares approximation, we assume that we have a function, g — g[x) defined over a finite set of 
points and we seek an approximate value of t; at one of these points, Xo, based upon the values at neighboring points. 
We assume that the function we seek is smooth enough to be expanded about the point Xq in a finite polynomial basis 

f x"^ x"" 1 

Thus, we can write 

rip 

f{x)^YaiPi{x-Xo) + f{xo) (17) 

where is a vector of coefficients. For a simple polynomial as above, the derivatives of / at Xo are then just the 
coefficients themselves. If we have a data point at Xo and we want the approximation to pass though the remaining 
data points, we can write a set of least-squares equations 
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= V'aiPi(xj - Xo), (18) 



where fj — f{xj) — f(xo), Pi(xj — Xq) is the ith polynomial basis member evaluated at Xj — Xo, and are the 
expansion coefficients. In other words, we seek a vector a which satisfies the linear equation: 

f = P-a (19) 

where /, P, and a are the terms in Eq. |l^ written in matrix/vector form. 

In order to solve for a, we need at least as many data points in the neighborhood about Xo as we have terms in the 
polynomial expansion. Rather than trying to pick out exactly the right number of points, the most efficient procedure 
is to simply select more data points than basis functions. Since the system of equations is overdetermined, P, is a 
rectangular matrix and we use must either singular value decomposition or other pseudo-inversion method to invert 
P. Furthmore, we can improve the stability of the least squares procedure by weighting each data point so that points 
farther away from the central point receive the least weight. Thus, for the "weighted" least squares procedure we 
solve 

j-n ^pn ,^ (20) 

where fp = {fj — f{xoj)ujj, Pfj = i^jPij, a is a vector of undetermined coefficients, and ujj is the weight assigned to 
point Xj. 

We have found that the convergence of the least-squares procedure is greatly improved by using a logarithmic form 
of the densityBEZN, 

P = e^(-), (21) 

where g is a polynomial of order Up. This representation is a useful way to transform a non- linear model to a linear 
one. For example, in computing the quantum potential we can write 

J-y^^^e-^/'v\<'/^ (22) 
VP 

Working through the derivative terms yields: 

,-,/2^2^,/2 ^ 1 + 1^2^^ (23) 

where the "comma" delimited subscripts in g.„ = dfj_g denote taking derivative with respect to coordinate fi. 

Finally we write the equations of motion for the velocity in terms of the derivatives of the various fields in component 
form: 

( 1 



ma^(x) = -a^ \V{x) ~ — (^g,,, + j (24) 

p{x, t + 5t) = e'^'^^'^'-'^^^'pix, t). (25) 

Some final notes regarding our implepientation of the MWLS scheme are also in order. We propagate trajectories 
using the sympletic Verlet algorithm cil. Convergence with respect to time step is checked by comparing the results 
of two 1/2 time steps to one full time step. If the differences between position, energy, velocity, or p computed via 
the two separate procedures varied by more than 1 : 10^, the shorter time step results were taken, 6t was reduced 
by a factor of 0.75 and propagation continued. However, if the longer time step results were sufficiently accurate, 5t 
was increased by a factor of 2 for the next iteration cycle. For the least squares fitting procedures, we found that 
a Hermite polynomial basis provided the most robust basis as measured by the for each fit. We also weighted 
the fits using a Gaussian weighting function centered about Xo adjusted so that the point farthest from Xo received 
a weight of 0.01. For the one dimensional results presented here, our code reproduces results obtained via Wyatt's 
MWLS particle code described in Ref.El 
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IV. EXAMPLE CALCULATIONS 



A. Harmonic Oscillator 



As a primary example, we consider the quantum trajectories for a harmonic osciUator and compare these to the 
classical trajectories starting from identical starting points. In Fig. 1 we show the time evolution of an ensemble of 
quantum trajectories for a harmonic system with period of r = 2tt/u} ~ 888. 57a. u. Unless otherwise noted explicitly, 
we shall use atomic units from here on, setting H = 1. The first two cases (A & B) are the quantum trajectories for a 
particle with mass m — 2000a. u and with mass m — 200a. u. respectively, each starting with identical initial quantum 
densities 

p{x) = e-'^t"-"")' (26) 

with Xo = 3.0 and /3 = 0.3. 

In each of the cases shown here, we used 100 points spaced evenly about the center of the initial density. In fourth 
plot (D) we show the trajectories obtain from a purely classical calculation involving the same 100 points. The effect 
of the quantum potential is clearly seen as the trajectories traverse the minimum of the well at a; = 0. For the heavier 
mass (in A), the trajectories tend to focus at the well minimum, just as in the classical calculation (in D). However, 
the trajectories do not cross each other as do the classical trajectories. This is a crucial characteristic of Bohmian 
trajectories. In essence, the density must be single valued and the Jacobian of the volume element carried by each 
point must remain positive definite for all times. Hence, if any pair of trajectories were to cross, these conditions 
would be violated. For case B, in which m — 200a. m , the quantum potential is lOx stronger than in case A and 
hence the trajectories become more diffuse as the particle traverses the bottom of the well. Interestingly enough, if 
we instead choose /? = 1/muj, the trajectory lines remain evenly spaced throughout the calculation as shown in C. 

This last result can be understood from the fact that the wavepacket used in Fig. l.C is a coherent state and hence 
the time evolution of the density is given by 

pix, t) = e-""^(^-^'' cos{ut)f ^27) 

where Xo is the original centroid of the wavepacket. The auantum potential can be easily derived from Eq. 26 and 
the action along any given Bohmian trajectory is given byEI 



(28) 



S{x, i) — —hcot — ^xxo sm{ujt) — —x^ sin(2cjt)^ . 

Thus, each trajectory is a solution of 

mx = dxS = —mujXoSiiT^ujt. (29) 

and belongs to the family of cosine curves of period lo: 

x,{t) = x^{Q) + Xo{cos{ujt) - 1). (30) 

I.e. the Bohmian particle undergoes simple harmonic motion of amplitude Xq about the point Xi{fi) — Xo. Furthermore, 
the energy along a given Bohmian trajectory is computed by taking the time derivative of the action (Eq. |2^ ): 

E{x, t) = dtS - ^muj'^{2xXo cos(tcj) - xl cos{2tLo)). (31) 

Note that this is not constant in time as would be expected for a purely classical system. The quantum force thus 
acts as a time dependent driving force on the ensemble of trajectories. If we were to have chosen the initial center of 
the density to be at exactly x — 0, then E = huj/2 and the quantum force would exactly counter-balance the classical 
potential force. In essence, the particles would remain fixed in their spatial positions. 

In the case of the quantum harmonic oscillator, the quantum trajectories computed using the MWLS procedure 
agree with the exact analytical results one can obtain for this system. Furthermore, the procedure appears to be 
stable for long periods of time. We have carred out calculations up to 10 oscillation periods without considerable 
accumulation of error. We now move on to the more formidable problem of tunneling in a double well potential. 
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B. Quantum Tunneling in a Double Well Potential 



Tunneling and barrier penetration represent phenomena which are unique to quantum mechanics. It is well under- 
stood from classical mechanics that if a particle lacks sufficient energy to surmount a potential barrier, the particle 
will be reflected and the barrier is impenetrable. However in quantum mechanics, a particle can "pass through" a 
barrier. This can be rationalized in a variety of ways. Perhaps the most consistent is that if the particle is in brought 
to a halt at the classical turning point, both its position and momentum are defined to exact precision, which would 
violate the uncertainty principle. In the quantum trajectory approach, the quantum potential acts to enforce this 
behavior by effectively "lowering" the barrier to permit trajectories to penetrate through the barrier and raises the 
barrier once sufficient amplitude has passed. To see how this comes about in the context of a particle based approach, 
let us consider tunneling through the double well potential: 

V{x) = ax'^ - bx^ (32) 

In this example, we take a = 0.007 and b = 0.01 in atomic units, which gives two tunneling states below the barrier 
sX E = — 369.827cm~^ and -313.918 cm~^ and a barrier height of Vf, = 786.24cm~^. The initial Gaussian density 
was centered about the center of the right side well at Xo = yjbjla with /? = \J 4bm corresponding to a harmonic 
approximation to a state localized in the right hand well. In the lower frame of Fig. 2 we show quantum trajectories 
for the initial crossing of the wavefunction from the right to the left side of the well and in Fig. 3 we show snapshots 
of the quantum density, the quantum potential, V + Q, and V{x) for various times. The trajectories cross from +x to 
—X corresponding to top to bottom in Fig 2. The initial expansion of the paths is primarily due to the expansion of 
the wavefunction in the well. Since the initial state is not in a stationary state of the potential, there is a net force on 
the particles due to the difference between the Q' and V. Once the particles have passed through the barrier region 
(at X = 0), they encounter the repulsive potential wall and are reflected. Note that the particles penetrate fairly 
deeply into the repulsive region. However, they carry very little net density. Unfortunately as far as our calculations 
are concerned, we run into numerical difficulties in computing the quantum potential in this region. As a result, the 
congruency condition fails to hold and trajectories begin to cross (after t = 650) and our results are meaningless 
beyond this point. 

However, prior to this point in time, we can compute the effective barrier height due to the quantum potential at 
a; = 0. In Fig. 4 we plot the effective barrier height, 14// = Q{0) -|- V;,, as a function of time. As indicated above, the 
quantum potential from the initial Gaussian distribution effectively "lowers" the barrier allowing trajectories to pass 
from the right to the left. However, once the initial expansion of the density is accomplished, the barrier begins to 
rise and the effective barrier height appears to approach an asymptotic value. 

If we assume that the lowest two eigenstates of the double well can be approximated as symmetric and antisymmetric 
linear combinations of ground state harmonic oscillators centered about the left and right minima in the well, we can 
compute the approximate quantum potential at a; = quite easily. Let us write the approximate wavefunction as 
(neglecting a common phase factor), 

^{x,t) = -^(e--V+(a;) + e+'^'^-{x)), (33) 

where (/>±(a;) are the symmetric and anti-symmetric tunneling states split by Hlo = {E-E+)/2. Taking, 

Mx)^^{Mx)±Mx)), (34) 

and using the approximation 

M^) = {^y\—^^''^^^^'/'\ (35) 
the quantum potential at the barrier is given by 

Q{0) = + \mu;lxl{cos{Uu) - 3)) (36) 

Thus, the effective barrier separating the right and left sides of the system is at its minimum when p is localized in 

the right or left hand well and greatest when pdelocalizcd in a 50:50 mixture of right and left hand states. Moreover, as 
the separation between the respective wells increases, the harmonic oscillator approximation to the lowest eigenstates 
becomes better and better. Correspondingly, the quantum potential becomes the parabola. 
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Qix)^^^^mu;l{x±Xof. (37) 

Since the quantum force and the classical potential force will be nearly equal and opposite, the particles will remain 
almost motionless in the initial well. 



V. PILOT WAVE TRAJECTORIES 



As noted above, the trajectories themselves can be obtained by propagating a solution to the time-dependent 
Schrodinger equation. In this "pilot wave" scheme, the quantum potential and equations of motion for the Bohmian 
particles may be computed directly from the wave function itself. Alternatively, if we have an accurate representation 
of V' and can compute v[ip] via Eq. |l|, we can obtain what should be a set of "exact" Bohmian trajectories associated 
with the quantum wavefunction-. To compute such trajectories for the double well example shown above, we used a 
discrete variable representationEa (DVR) of Gauss-Tchebychev quadrature points to represent both the wave packet 
and the time-evolution operator on a finite spatial grid of 200 quadrature points spanning 2.5 Bohr on either side of 
the barrier. For dynamics at low total energy, this representation is certainly more than adequate as we were able to 
converge the lowest 20 eigenstates of this well to 1 : 10^. Because the instantaneous position of a Bohmian particle 
at x{t) will not generally correspond to the position of a quadrature point, the velocity can not be computed directly 
from the DVR wavefunction, which is only defined on the quadrature points them selves. Rather, we transform ip 
from the DVR to a finite basis representation (FBR) using the unitary transformation, 

^FBB. ^ ^ . ^DVR^ (38) 

where U is the unitary transformation associated with a A^-point Gauss-Tchebychev DVR. 



The elements of the vector tp^^^ are the expansion coefficients of the wavefunction in a finite polynomial basis. 

n 

V'(x,t)=^^™(OT,(x). (40) 

i=l 

For the Gauss-Tchebychev quadrature scheme used here, basis functions are: 

Ti{x) = y^sin(i7r(a; - A/2)/A). (41) 

Eq. 37 gives the wavefunction at any point x G [—A/2, A/2] and we can use this information to compute the 
velocity of a particle guided by ip. In the top frame of Fig. 2 we show the Bohm trajectories for the double well 
system choosing identical conditions as above, this time computed via the DVR pilot wave based procedure. For 
clarity we show only those trajectories for particles starting to the left of the initial centroid of the density. At long 
times, the trajectories transmitted to the left hand well tend to bunch together. This is due to the fact that the finite 
basis functions in Eq. 38 go to zero at a; = ±A/2 and force the system to have an artificial node at the boundary 
points. Consequently, the quantum potential is infinitely repulsive at the end of the grid and all particle trajectories 
are repelled from these regions. For practical purposes, the density carried by these trajectories is negligible and the 
overall effect on the wavepacket dynamics is minimal. 

Notice, however, that between t — 600 and t = 800, a number of trajectories near the barrier at a; = are 
deflected from each other. This is due to the fact that p develops oscillatory structure due to constructive and 
destructive interferences between on coming and reflected components of the wavefunction near the barrier. The 
density eventually develops a node at about t = 900 a.u. Surprisingly, and disturbingly, this behavior is absent in 
the MWLS trajectories. In Figure 5 and Figure 6 we compare the evolution of p using the DVR (p^^^) and MWLS 
^pMWLS^ methods. Initially, the agreement is quite good. At longer times the agreement is very poor with p^^^^ 
even "looping back" through itself. At intermediate times, as p is reflected by the repulsive barrier in the —x half 
of the double well, oscillations develop in p^^^ whereas no such oscillations are seen in p'^^^S ^ gijj(jg ^jjg quantum 
potential is a measure of the local curvature of p, these oscillations translate into attractive and repulsive regions 
in the quantum potential. Eventually as p^^^ goes on to form nodes (at approximately t = 700 au and 900 au 
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trajectories are guided away from such regions via the quantum force. This accounts for the significant deflections in 
the DVR trajectories but does not tell us why this fails to occur in the MWLS case. 

The answer is basically that of sampling resolution. As the MWLS trajectories evolve, they by in large tend to 
spread apart. This restricts our ability to resolve any structure in p^^^^^ which occurs on lengthscales finer than 
the distance of separation of the MWLS trajectories at a given point in time. Hence, the oscillatory structure seen in 
the DVR calculations (e.g. at ^=550 au) are too fine to be resolved by the MWLS trajectories. Thus, the quantum 
potential will be too smooth in this region. The DVR calculations are also band width limited due to the finite 
spacing of the DVR points. However, in the DVR grid used in these examples, this spacing is fixed at 5x = 0.049 
Bohr, where as in the MWLS case the spacing between trajectories in the region of the barrier ai t = 550 is w 0.10 
Bohr, which slightly courser than the structure seen in p^^^ in this region. Consequently, the MWLS trajectories 
lose band width resolution precisely where it is needed the most in this case. 

Let us compare the DVR and MWLS based trajectories on a head to head basis. Four sets of trajectories extracted 
from Fig. 2 are shown in Figure 7. First, for trajectories originating near the leading edge of the MWLS grid, we 
expect that the two results should deviate rather early due to the fact that the fitting procedure used to compute the 
quantum potential using the MWLS is least accurate near the end of the grid (Trajectory # 9). Surprisingly, however, 
the agreement is quite good. In fact, it is appearent that the deviatiation is more likely due to the artificial boundary 
conditions imposed by the basis function used to represent the quantum wave packet in the DVR calculations. For 
trajectories originating near the maximum of the initial wavepacket (Trajectory #50), we see quantative agreement 
between DVR and MWLS based trajectories since p is smooth function of x in this region and remains so over the 
course of both calculations. However, let us compare the two trajectories labeled #38 and #39 from each calculation. 
These two trajectories originate side by side separated by 6x = 0.02Bohr. In the DVR calculation, these trajectories 
are defiected by the oscillations in p with #38 deflected slightly towards the —x direction and #39 in the +x direction. 
These deflections are not observed in the corresponding MWLS trajectories. 



VI. DISCUSSION 



In this paper, we present our implementation and computational analysis of the de Broglic-Bohm hydrodynamic 
equations using a particle based description. The formalism itself allows an elegant interpretation of quantum dynamics 
in geometric terms. Instead of wave functions, we have geometric ray lines. Surprisingly, this intuitive connection 
between quantum wave mechanics and geometry has not made it into the "main stream" even though the original 
seeds for this interpretation can be found at the very advent of quantum theory. However, we do note that this 
view point is seeing a resurgence in popularity as measured by the number of papers which have used the Bohmian 
construction in vastly different areas of quantum physics. 

What is important, however, to point out that we are only showing the best case scenarios in this paper. We have 
found that trajectories tend to cross cases in which p is very small, trajectory lines are far apart, or when the potential 
is significantly anharmonic. Some of this is exhibited in Figures 2 and 3. This is most certainly a result of numerical 
instabilities in the MWLS procedure which we have elaborated upon in this paper. 

We also run into particular dificulties in dealing with nodal points in the density, in other words, when p{x) — 
'ijj{x)\'^ = 0. At such points, it is not possible to fit the coefficients of a polynomial expansion using points on either side 
of the node since the function does not have continuous derivatives at the node. We have tried various computational 
alternatives to dealing with such points. One alternative is to use only trajectories and densities on one side of the 
node or the other when computing Q{x). Under this approach we compute derivatives by approaching the node 
from the right and from the left. Another alternative we have tried is to introduce a log(x — x„) basis function 
for stars which include nodal points. Both of these approaches give very good results for cases in which the initial 
wavefunction was taken as an odd parity eigenstate in a symmetric well. However, both cases require prior knowledge 
of the location of the node and that parity be a constant of the motion so that the node does not move. For cases in 
which parity is not a constant of the motion, nodal points can occur over the natural course of the evolution of the 
quantum wavefunction, even if the initial wavefunction has no nodes initially. We speculate that a solution may be 
to construct 

R^ = Vple'^" (42) 

along each ^^(i) JLrajectory where p — 0, 1,2, . . . accounts for the branch cut taken when evaluating on either 
side of a node. c3 In essence, on one side of a node R > and p ~ 0, 2, while on the other side, R < and 
p = 1,3, .. . This would preserve analyticity across the node allowing us to calculate Q using a simple polynomial 
basis.i-^jVt present, the treatment of nodes and their evolution is an open problem which we will take up in a future 
work.cZl 
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FIG. 1. Quantum trajectories for harmonic oscillators. In each case,the oscillation period, r = 27r/a' = 888.57oM. In cases A 
& D ,m = 2OOO0M while in case B, m = 200om. Case D is a set of classical trajectories {Q = 0) for this system. 
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FIG. 2. Bohmian trajectories computed via wavepacket propagation (DVR) versus particle based trajectory method 
(MWLS). 



10 



X 



DVR 





FIG. 3. Evolution of quantum density, p (dots), quantum potential (short dashed lines), V{x) (solid line), and 
V + Q{longdashedlines) for a proton in a double well potential at various times. 
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FIG. 4. Effective barrier relative to bottom of well encountered by quantum trajectories in the double well problem in 
crossing from right to left hand well. 
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FIG. 5. Evolution of the quantum density, p, as computed via MWLS based methods. 
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FIG. 6. Evolution of the quantum density, p, as computed via DVR based methods. Note in particular the interference 
structure which develops in the DVR calculation which does not appear in the MWLS trajectory calculations shown in Figure 
5. 
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FIG. 7. Comparison of 4 trajectories computed via DVR (dashed) versus MWLS (solid) methods. 
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